Heat transport in model jammed solids 
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We calculate numerically the normal modes of vibrations in 3D jammed packings of soft spheres as 
a function of the packing fraction and obtain the energy diffusivity, a spectral measure of transport 
that controls sound propagation and thermal conductivity. The crossover frequency between weak 
and strong phonon scattering is controlled by the coordination and shifts to zero as the system is 
decompressed towards the critical packing fraction at which rigidity is lost. Below the crossover, 
the diffusivity displays a power-law divergence with inverse frequency, which suggests that the 
vibrational modes are primarily transverse waves, weakly scattered by disorder. Above it, a large 
number of modes appear whose diffusivity plateaus at a nearly constant value independent of the 
inter-particle potential, before dropping to zero above the Anderson localization frequency. The 
thermal conductivity of a marginally jammed solid just above the rigidity threshold is calculated 
and related to the one measured experimentally at room temperature for most glasses. 

PACS numbers: 45.70.-n, 61.43.Fs, 65.60.-(-a, 83.80.Fg 
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I. INTRODUCTION 

The thermal and mechanical properties of disordered 
solids can differ dramatically from those of crystalline 
materials Prominent among the anomalous proper- 
ties are the specific heat and thermal conductivity, which 
display common features at sub-helium to room temper- 
atures, in amorphous materials ranging from glasses to 
plastics and even frozen grease (2|. This commonality 
suggests that the explanation of the unusual features of 
disordered solids may involve general physical principles 
that transcend detailed information about the chemical 
structure of specific compounds [1, |3, @ • 

This paper focuses on the intermediate temperature 
regime, IK < T < Tioom- In this regime, the ther- 
mal conductivity has a plateau followed by a nearly lin- 
ear rise at higher T, in contrast to the sharp drop seen 
in crystalline materials in the same range [6|. In the 
temperature range of this thermal-conductivity plateau, 
the ratio of the heat capacity C(T) to the expected T^- 
dependence predicted by the Debye model for crystalline 
solids, C{T)/T^, exhibits a prominent peak, termed the 
boson peak, which is a hallmark of amorphous solids 
0,0. 

At lower temperatures, T < IK, the thermal conduc- 
tivity exhibits a T^-rise, in contrast to the dependence 
observed in crystals Q. Meanwhile the heat capacity 
increases linearly in T in contrast to the rise 
predicted by Debye for long- wavelength sound modes. 
These low-temperature features have been rationalized 
by invoking the scattering of long- wavelength phonons off 
two-level systems, posited to arise from groups of atoms 
tunneling between two minima 0, • While this view 
of the low-temperature regime is widely accepted, little 
consensus exists regarding the intermediate temperature 
regime. In particular, the origin of the plateau in the 
thermal conductivity is controversial, and the question of 
whether the plateau in the thermal conductivity is linked 



to the boson peak remains unresolved. 

In this article we present a systematic study of the pres- 
sure (or packing fraction) dependence of the energy diffu- 
sivity of random packings of soft frictionless spheres. The 
energy diffusivity, d{uj), quantifies how far a wavepacket 
narrowly peaked at a frequency u can propagate. We 
study a class of models in which spheres interact via a 
repulsive potential that vanishes at a well-defined dis- 
tance that defines their diameter. These models possess 
a j amming / uni amming transition as a function of pack- 
ing density Thus, the strength of the elastic moduli 
can be tuned continuously downwards by decreasing the 
density of the packing until the particles no longer inter- 
act. At that critical density, the disordered solid unjams 
to form a liquid. As a result, the onset of the excess vi- 
brational modes uj* can be pushed to zero by decreasing 
the packing fraction [l^. Numerical calculation of the 
density of states and the diffusivity enables us to track 
the pressure dependences of the boson peak frequency 
UJ* and the transport crossover frequency ujd individually 
and to compare them. 

We find that the diffusivity [13] displays a well-defined 
kink at a frequency Ud that separates the low uj regime 
of divergent diffusivity from a characteristic plateau that 
persists up to high u where Anderson localization sets in. 
Moreover, ujd and the boson peak frequency, to*, are not 
only comparable in magnitude, but decrease in tandem 
as we decrease packing fraction towards the unj amming 
transition. This shows that the transport crossover is 
linked to the excess modes. 

A unique feature of disordered sphere packings is that 
the transport crossover can be understood perturbatively 
by means of simple scaling laws parameterized in terms 
of the distance from the critical unjamming transition, 
denoted as Point J [llj. Furthermore, some of the dis- 
tinctive properties of amorphous solids are manifested 
in their most extreme form as the unjamming transition 
is approached. For example, in that limit the diffusiv- 
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ity plateau extends all the way to the lowest frequencies 
studied in the simulations. This suggests that the origin 
of the diffusivity plateau can be traced to properties of 
the unjamming transition p^ . 

The outline of this paper is as follows. In Section 
im we provide the necessary background on the link be- 
tween vibrational dynamics and heat transport on which 
our work is built and review the scaling properties of 
jammed solids. In Section lllli we review the methodol- 
ogy adopted to calculate the energy diffusivity using the 
Kubo formula, which enables a first-principles calcula- 
tion of the diffusivity for computer-generated packings. 
In section llVi we present our results for the transport 
crossover and show that this occurs at the boson peak 
frequency. We also present a scaling analysis that ra- 
tionalizes how the main features of the diffusivity de- 
pend on the distance from the unjamming transition for 
both Hertzian and harmonic interactions. Section |V] fo- 
cuses on the most striking transport signature of jammed 
packings: a plateau in the diffusivity as a function of fre- 
quency. Such a plateau was postulated by Kittel 
in order to explain the temperature dependence of the 
thermal conductivity of amorphous solids at intermedi- 
ate temperatures; our results provide evidence that this 
scenario is correct. We explain the origin of the diffusiv- 
ity plateau, starting from a set of assumptions concerning 
the nature of the vibrational modes above the transport 
crossover. The AC thermal conductivity at point J is 
then obtained in section IVTl In Section IVlII wc conclude 
by summarizing the broader message of this article: the 
thermal conductivity of various amorphous materials un- 
der pressure can be explained from the vibrational modes 
at Point J, which controls energy transport at higher den- 
sities in the manner expected for a critical point. 



II. BACKGROUND 

It was originally suggested by Kittel [l3| that the in- 
termediate temperature properties might be understood 
from the microscopic vibrational dynamics of amorphous 
materials. For instance, the boson peak at T* implies 
the presence of vibrational modes in the density of states 
at oj* = ksT* /h in excess of the usual Debye counting. 
Similarly, the plateau in the thermal conductivity could 
result from a crossover in the behavior of the mean free 
path of phonons [3 ■ 

Several theoretical models have been advanced to ex- 
plain the boson peak, starting from distinct physical 
mechanisms such as the existence of resonant (quasi)- 
localized modes [Tol . [20l . [2ll [2^ , anharmonic interactions 
induced by the presence of defects [13, [1^, the breakdown 
of continuum elasticity below a characteristic length scale 
15, [2g| or quenched disorder in the elastic constants 
27, l28l [29|. Some of these models find that the onset of 
the excess vibrational modes coincides with a crossover 
from weakly-scattered plane waves to strongly-scattered 
vibrational modes that are delocalized and poorly con- 



ducting [so, |3l|, suggesting a connection between the 
boson peak and the plateau in the thermal conductiv- 
ity. 

The connection between the boson peak and the trans- 
port crossover has been probed using Brillouin scat- 
tering measurements but no firm conclusion has been 
reached to date [H, [s^. A recent experimental study 
suggested that the excess modes in the density of states 
appear around the loffe-Regel frequency 3j|, at which 
the mean free path of longitudinal phonons approaches 
their wavelength [2^. This heuristic criterion to esti- 
mate the crossover between weak and strong scattering 
was also used in an independent study that challenges 
the previous claim by concluding instead that the loffe- 
Regel limit takes place at a fre quen cy significantly higher 
than the Boson peak frequency [35| . Classic experimental 
studies of Raman scattering point to another important 
piece of experimental evidence, namely that the vibra- 
tional modes at the Boson peak are transverse in char- 
acter [sH, [131 ■ This observation is supported by simula- 
tions of silica [ssl, '39'| and soft sphere glasses [ioj . Recent 
numerical simulations provided additional evidence that 
suggests the equality of the Boson peak frequency with 
the loffe-Regel limit for transverse phonons • 



A. Vibrational dynamics and heat transport 

Instead of focusing directly on the heat capacity, C(T), 
and thermal conductivity, n{T), we consider the density 
of states, -D(w), and diffusivity, A heuristic argu- 

ment for the relation between d{ijj) and k{T) is as follows. 
For a system in a temperature gradient, it is well known 
that the heat diffusivity obeys the relation d — kV jC . 
Thus, K — dC /V . This relation can be generalized mode 
by mode. Thus, the heat capacity is [4 ^S,, ,41j 



C(r)=^C(a;,,r) 



(1) 



where the sum runs over all vibrational modes i and 
C(w,;, T) is the heat capacity per mode, which is obtained 
from the Bose-Einstein distribution and is a universal 
function that characterizes the heat carried by a mode 
of frequency uji at temperature T. Similarly, the thermal 
conductivity is 



(2) 



where d{<jJi) is the energy diffusivity of mode i. 

We may recast Eqs. [l][2] in continuum form using the 
density of vibrational states: 



C{T) 



AujD{uj)C(uj,T) , 



1 r 



du)D{u))d{u))C{uj,T) . 



(3) 
(4) 
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Both D(u;) and d(u}) are strongly structure- dependent: 
the density of states and diffusivity are the fingerprints 
of the vibrational modes in the solid and control its heat 
capacity and thermal conductivity. 

Inspection of Eq. ([3]) reveals that the prominent boson 
peak observed at the temperature T* in most amorphous 
solids is triggered by a large number of excess vibrational 
modes that show up in the density of states at a char- 
acteristic frequency uj* ^ It is known empirically 
that ui* increases as the sample is compressed [42,li3|, a 
property shared by the soft sphere packings investigated 
in this study. By analogy, Eq. (jH) suggests that the ori- 
gin of the thermal conductivity plateau around T* can be 
similarly traced to the existence of a transport crossover 
in d(uj) at uj*. 

At very low frequency, the diffusivity can be factored 
out as the product of the speed of sound v and £{u}) [1]: 

d{co) = ^e{u). (5) 

The mean free path £{uj) typically diverges as cj — > 
because the corresponding vibrational modes can be re- 
garded as long wavelength plane waves weakly scattered 
by disorder. However, once the mean free path shrinks to 
the order of the wavelength, Eq. ([5]) breaks down j23l.l43|. 



B. Jammed sphere packings 

We study a model of amorphous packings of frictionless 
spheres interacting via the repulsive pair potential V{rij) 

V{rij) = if > (7ij , (6) 

where the distance between the centers of particles i and 
j is denoted by and the sum of their radii by cTy . 
We generate T = packings by conjugate-gradient en- 
ergy minimization according to the procedure described 
in Ref. [ll|- Irrespective of the value of a, this model sys- 
tem exhibits a jamming/ unjamming transition at T = 
at a packing fraction cj) — (pc (Point J) at which the par- 
ticles are just touching each other and there is no over- 
lap "r^. 

The zero-temperature jamming/unjamming transition 
has mixed character. At this transition, the average coor- 
dination number, z, jumps from zero to the min- 
imum value required for mechanical stability, the "iso- 
static" value Zc = 2D [4^ , where D is the dimensionality 
of the sample. At densities lower than (pc, particles are 
free to rearrange while above (j)c at = (j)—(j)c , the sys- 
tem behaves as a weakly-connected amorphous solid with 
an average coordination number that scales as a power 
law with an exponent consistent with 1/2 [Til, lisj: 

Az = z~ . (7) 



In addition, both elastic moduli exhibit scaling behavior 
near the jamming point consistent with [Til. IZq: 

G -A0"-3/2^ 

B -A(/)"-2 . (9) 

For harmonic repulsions (a = 2), the bulk modulus is in- 
dependent of compression while the shear modulus van- 
ishes as point J is approached. The bulk modulus scales 
as the second derivative of the potential with respect to 
compression, while the scaling of the shear modulus does 
not follow this naive scaling. 

The ratio G/B ~ Acffi-^ of the two elastic moduli is in- 
dependent of a and controls the relative contribution of 
transverse to longitudinal waves at low frequency. This 
can be checked explicitly by considering that the phonon 
density of states at very low frequency satisfies the ubiq- 
uitous Debye law D{uj) ~ ^ except at Point J where, 
as we shall see, the Deybe regime is completely swamped 
by a new class of vibrational modes. The transverse and 
longitudinal speeds of sound vt and vi are proportional to 
the square root of the shear and bulk moduli, respectively 

Vt ^Acj)^, (10) 
VI ^A(p^ . (11) 

Upon substituting into the Deybe formula, D{uj) ~ 
Lu'^/v^, Equations pU|) and (fTT|) imply that the ratio of 
the transverse to the longitudinal density of states at low 
frequency, Dt/Di ^ A(j)~^^'^, becomes arbitrarily large 
as A(j) 0. Thus, the density of states is dominated 
by transverse modes at low frequencies where wave-like 
behavior is expected, irrespective of the potential. 

Numerical studies [Til [l2| have revealed the presence 
of excess vibrational modes that contribute to a plateau 
in the density of vibrational states above a characteristic 
frequency, uj*. Close to the jamming point u* increases 
with density consistent with the power law [l^ 

UJ* - Ac/}"^ . (12) 

Thus, the plateau extends to zero frequency in a 
marginally-jammed solid (i.e. a packing of particles just 
above the onset of mechanical rigidity) . 

III. METHODS AND MODEL 

The energy diffusivity, d{uj), introduced in Sec. Ill Al 
can be viewed physically in terms of the behavior of a 
spherically-symmetric wavepacket narrowly peaked at a 
frequency lu and localized at position r at time t — 0. 
The spatial width of the wavepacket spreads out because 
the normal mode components from which the wavepacket 
is constructed propagate in all directions. The time- 
independent diffusivity, d{uj), is given by the square of 
the width of the wavepacket at time t, divided by t at 
long times t 
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If the width grows linearly in time, the difFusivity is in- 
finite; this corresponds to ballistic propagation. If width 
grows with the square-root of time, a finite difFusivity is 
obtained; this corresponds to diffusive propagation. As 
first noted by Anderson, a third possibility exists, namely 
that the width of the wave-packet saturates to a constant 
value over which the vibration is localized [47l|. Such lo- 
calized modes, typically occur at high to [4l|. The diffu- 
sivity is vanishingly small and d{ijj) cannot be factorized 
into the product of and v, as in Eq. [51 because no 
speed of sound can be associated with such vibrational 
modes. 

In this study, we calculate the diffusivity by evaluat- 
ing the Kubo formula for d{uo) directly [3] for computer- 
generated packings in terms of the normal modes over 
the entire frequency range available. The rationale be- 
hind this choice is two-fold. On one hand, we use the 
energy diffusivity as a spectral measure of transport to 
probe the character of the vibrational normal modes of 
a jammed solid. On the other hand, we use the jammed 
solid as a model amorphous structure whose transport 
properties can be studied as a function of pressure sim- 
ply by varying the density relative to that of the un- 
jamming transition. This allows variation over orders of 
magnitude of pressure, which cannot be realized in more 
realistic models of molecular or network glasses. 



A. Review of the Kubo formula for the energy 
diffusivity 

The analogy between the thermal conductivity in the 
phononic problem and its better-studied electrical coun- 
terpart underlies many of the mathematical techniques 
and physical concepts employed in our investigation. 

The Kubo formula for the energy diffusivity was de- 
rived in a convenient form by Allen and Feldman in Ref. 
0. Consider the volume — averaged energy current S 
that arises in response to an applied thermal gradient 
VT. In linear response, S is given by 



S = -kVT, 



(13) 



h that obeys the continuity equation: 

^ = -v..(^), 



(15) 



where s{r) is the local energy- flux operator (assumed to 
be isotropic for simplicity) prior to the spatial average 
that leads to S. In the case of the Einstein relation, the 
conserved quantity is the particle number and the role of 
s{r) is played by the particle current.) 

Our aim is to extract an expression for the diffusivity 
d{uj) by comparing Eq. ([14]) to Eq. (g]) when the DC 
limit of K(r, il) is taken, that is, when ^ 0. We first 
recast Eq. (|14p in terms of a discrete sum over modes i: 



{Sfj,)ij{S^)jiS{uji-ujj-Q) , 



h{u!i — UJj) 



(16) 

where rii is the equilibrium occupation number for bosons 
rii — [exp{/3huji) — 1]^^ and («S'^)ij is the matrix element 
of the energy flux operator in the /x direction. This ma- 
trix element can be computed from the vibrational nor- 
mal modes, which are obtained from the dynamical ma- 
trix , -ff™, whose i*'' normalized eigenvector is denoted 
by ei{m] a), where {m, n\ and {a, /?} label particles and 
their Cartesian coordinates respectively [48]. For disor- 
dered solids, the modes must be determined by numerical 
diagonalization of the dynamic matrix. 

In an isotropic system a scalar thermal conductivity 
K,{T) can be meaningfully defined from the trace of the 
tensor k^u{T) 



(17) 



To simplify Eq. (jl6p further, consider that in the limit 
^ the delta function forces the factor {nj — ni)/(u)i — 

uoj) to become —dru/duJi. One can then use Eq. p7)) in 

conjunction with the identity 



VT 



where n is the thermal conductivity. More generally, the 
AC thermal conductivity k^u{T, J7), which relates the en- 
ergy flux in the /i direction to the time-varying tem- 
perature gradient in the v direction, di^Te^^^, is 

K^AT,n) = — dx dt e^(^^+''')* {s^{~ih\)s,{t)) , 

(14) 

where V is the volume of the system and the angular 
brackets denote an equilibrium average of the autocorre- 
lation function of the energy current operator S. 

(Note that Eq. is analogous to the Einstein rela- 
tion for the diffusion coefficient of an ensemble of Brown- 
ian particles in terms of the velocity-autocorrelation func- 
tion. The main difference lies in the fact that in the ther- 
mal problem the conserved quantity is the energy density 



(18) 

to show that the thermal conductivity of Eq. can 
indeed be factorized as indicated in Eq. ^ with d{uJi) 
given by [1] 

d{u;,) = ^Y.{nuj,)-^ \S,,\^5{oj,-u;j), (19) 
j 

where the matrix elements Sij read Q 



4 CUi LOj 



In the limit LOi — s- ujj enforced by the delta function in 



Eq. (fTO]) . the prefactor 



4 LJi UJj 



1. However, taking 



this limit requires special care when the Kubo formula 
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(derived in the continuum limit ) is evaluated for a fi- 
nite and isolated system with a discrete spectrum [4^. 
Inspection of Eq. ((20|l reveals that the diagonal matrix 
elements Su vanish. On the other hand any contribution 
to d{uji) coming from the non-diagonal matrix elements 
Sij with i =^ j is given zero weight when the delta func- 
tion S{uji — LUj) is strictly enforced. This difficulty can 
be circumvented by smoothing out the ^-functions in Eq. 
(fT9)l with the small finite width rj. 



TT[{u!i - UJj)^ + 



(21) 



This heuristic procedure is expected to give the correct 
"bulk" result as long as the broadening 77 is (a) much 
larger than the average level spacing, A, and (b) much 
smaller than any characteristic frequency scale relevant 
to the problem. In this paper, the broadening of the 
delta function rj in Eq. (pTjl is typically chosen to be 
approximately five times larger than the average level 
spacing, A. We have verified that our numerical results 
do not depend on this choice as long as conditions (a) and 
(b) are met. In the Landauer formulation of transport, 
it is not necessary to introduce the level broadening ij by 
hand because the inherent couplingof the system to the 
reservoirs plays an analogous role [50| . 

One advantage of studying the energy diffusivity in- 
stead of the thermal conductivity is that d(uj) is finite 
at nonzero frequency when evaluated at the harmonic 
level. By contrast, k{T) is generally infinite if anhar- 
monic corrections are ignored. This is because d(uj) in 
the integrand of Eq. [4] diverges too strongly at low lu 
due to phonons that are progressively less scattered with 
increasing wavelength. 

One main result of this paper will be that a marginally 
jammed solid is an exception to the rule that the thermal 
conductivity should diverge within the harmonic approx- 
imation. However, as the system is compressed above 
Point J, k{T) again diverges, as in the standard case. 

In order to cure this divergence, additional scattering 
mechanisms, beyond harmonic theory, are typically in- 
voked resulting in an additional contribution to the dif- 
fusivity, dc{Lu). Upon adding dc{uj) to the harmonic con- 
tribution, d{uj), (for example, as if they were two con- 
ductors in series HH), one obtains the total diffusivity 
dxitLi)' 



driuj)'^ = d{uj)~^ + dc{Lo)~ 



(22) 



The cut-off contribution ddtu) necessarily dominates at 
low temperature, but much progress can be made in un- 
derstanding the plateau in the thermal conductivity and 
the subsequent rise by studying the harmonic contribu- 
tion, which dominates as T increases. 



a 50/50 bidisperse mixture comprised of 250 < N < 
10,000 frictionless spheres with a diameter ratio of 1.4, 
interacting with potentials described in Eq.[n]with a = 2 
and a = 5/2. The packing fraction at the onset of jam- 
ming, (f>c, is characterized by the onset of a nonzero pres- 
sure. We determine 4>c and obtain T = configurations 
at controlled A(f> = — 0c as in Ref. 12]. For each con- 
figuration, we diagonalize its dynamical matrix and find 
the eigenvectors and the corresponding eigenfrequencies, 
which are measured in units of y'e/MCT^, where M is the 
particle mass (l^ . 

We primarily consider an "unstressed" version of this 
model which is particularly tractable. We use 

energy-minimized configurations obtained from numer- 
ically generated jammed sphere packings as described 
above. We then replace the interaction potential, V{rij), 
between each pair of overlapping particles with an un- 
stretched spring with the same stiffness, V"(rlJ), where 
is the equilibrium distance between particles i and j. 



eq 



Note that since 



eq 



takes a different value between dis- 



tinct particle pairs there will be a distribution in the local 
values of the elastic constants. On the other hand, the 
fact that all springs are unstretched guarantees that there 
are no net forces between particles in their equilibrium 
positions so that stable configurations for the stressed 
system are also stable in the unstressed one. 

The unstressed packings correspond to dropping terms 
depending on the first spatial derivative of the potential, 
V' , in the dynamical matrix obtained from expanding the 
energy around the equilibrium position of the particles 



6E = 



(23) 

where n, m index particles. 

The approximation of dropping the V stress term in 
Eq. [53] generates an interesting disordered system in its 
own right. The resulting off-lattice model, comprised of 
point particles interacting with relaxed springs, exhibits 
both spatial fluctuations in the local elastic stiffness as 
well as topological disorder (e.g., fluctuations in the local 
coordination number). Moreover, its amorphous struc- 
ture can be varied by changing the volume. 

The effect of the stress term can be seen from Eq. (|23l) . 
Since V is negative for repulsive interactions, the stress 
term lowers 6E and hence the mode frequency. We will 
discuss the effect of the stress term on the diffusivity of 
a jammed solid in Section IV CI 



IV. ENERGY TRANSPORT CROSSOVER IN 
MODEL JAMMED SOLIDS 



B. Model 

Our simulations are carried out on jammed sphere 
packings, as described in Sec. Ill B[ Speciflcally, we study 



Figure [T] shows a scatter plot of the mode diffusivity 
d{u}i) obtained from evaluating numerically Ea. p9ll2m) 
at four packing fractions At/) — 0.3 (green), 0.1 (red), 
0.05 (blue), 0.02 (yellow) for a 2000-particle packing. 
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FIG. 1: (Color Online) Plots of diffusivity, measured in units 
versus frequency, measured in units of e/a'^M for 
an unstressed packing of 2000 particles with harmonic inter- 
actions at packing fractions Af/i = 0.3 (green), 0.1 (red), 0.05 
(blue), 0.01 (yellow). The black dots indicate the crossover 
frequency ujd at each A0, while the dashed lines show a power- 
law of tJ"*, expected for weakly scattered plane waves. The 
plateau diffusivity is do ~ 0.35. The inset shows the packing 
fraction scaling of Ud- 



riodic system and the speed of sound v is the transverse 
one for most low oj modes near 0c j see Eq. (|10j|lip . In 
the continuum limit we expect that their density of states 
at very low a; will be given by the Deybe law Diu)) ~ 

The intermediate frequency regime is the one of most 
direct relevance to the intermediate temperature behav- 
ior of the thermal conductivity. Strikingly, this regime is 
characterized by a diffusivity, henceforth labeled as do, 
that is nearly independent of frequency (Fig. [T]). The 
notion of a frequency-dependent diffusivity has a long 
history dating back to Kittel's observation [l3| that the 
experimental curve for n(T) in many glasses at room 
temperature could be interpreted in terms of a nearly 
frequency-independent mean free path of the order of a 
molecular length. 

The onset of the plateau in the diffusivity is marked 
by a crossover frequency, that exhibits a peculiar scaling 
with the packing fraction A0. We henceforth label it as 
LOd- The inset of Fig. [1] reveals that Ud ~ A(/)°-^ for a sys- 
tem composed of harmonic springs. This is the same scal- 
ing with observed for the frequency u* above which 
a large excess of vibrational modes have been observed 
in previous studies of the density of states (see Eq. (fT^ 
with a = 2). 



Three distinct transport regimes can be clearly identi- 
fied in each curve, corresponding to ballistic, diffusive 
and localized vibrational modes. The crossover frequen- 
cies that mark the ballistic to diffusive crossover for each 
A(/), indicated in Fig. [1] as black dots, do not depend on 
N for systems of this size or larger. In what follows, the 
mode diffusivity data is presented as scatter plots of sin- 
gle particle configurations for clarity. We have tested that 
performing a frequency binning followed by a disorder- 
average over several distinct configurations confirms our 
conclusions p^ . 

At very high frequencies, Anderson localization sets 
in and the diffusivity drops rapidly as a result of local- 
ization of the vibrational modes. The contribution of 
localized modes to the thermal conductivity is negligi- 
ble if anharmonic effects such as hopping are ignored. 
Fig. [H shows that the localization frequency for particles 
interacting with a repulsive harmonic potential does not 
depend strongly A(/) close to the jamming point. 

At low frequencies, the energy diffusivity exhibits a 
strong frequency dependence characteristic of vibrational 
modes that are essentially phonons weakly scattered by 
disorder. As a comparison we have drawn black dashed 
lines in Fig. [T] indicating the power law divergence with 



A. Dimensional analysis 

Fig. [H shows that d{uj) is characterized by a well- 
defined crossover from ballistic to diffusive behavior. Our 
aim in this section is not to provide a rigorous derivation 
of the functional dependence of the diffusivity on fre- 
quency, but rather to understand how the defining fea- 
tures of the diffusivity (the height of the plateau, do, and 
the scaling of the crossover frequency, uJd) depend on ap- 
plied pressure. This is done by keeping track of how the 
fundamental parameters that enter the definition of d{uj) 
scale with the packing fraction A(^. First, however, we 
must understand how these parameters depend on the 
dimensional parameters in our system: the particle di- 
ameter, cr, the particle mass, m, and the energy scale for 
the potential, e. 

The diffusivity has dimensions of lirlMl^^ xhe natu- 

time 

ral unit of frequency in a vibrational system, by which the 
Lo axis in Fig. [l]is measured, is \J'^ where k is the hare 
elastic coupling of the solid. More generally an effective 
spring constant, fee//, can be defined by differentiating 
twice the interaction energy V{rij) of Eq. ([6]) and cval- 



expected for Rayleigh scattering of plane waves inci- uating the result at the average equilibrium bond length 



dent on uncorrelated scattering centers. Close inspection 
of the scatter plot reveals that the low lo peaks occur 
close to the discrete frequencies allowed in our cubic sim- 
ulation box of size L by the linear dispersion: 



< r- > of interacting neighbors: 



52l/(r„„) 



2t:v 



Qj,2 

^ nni 



(25) 



(24) 



where {p, q, r} denote the quantum numbers for the pe- 



There is a simple linear relation between the relative 
change of < rj^^ > upon compression and the corre- 
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spending macroscopic change of volume 



fj- < r^'? > 

^ ^ nrn ^ 



A<j) 



(26) 



where the numerical prefactor in the right hand side of 
Eq. [^was checked numerically [6^. 
From Eqs. (|25ll26p we see that 



a - 1 /A0 



(27) 



For harmonic repulsions, a = 2 and k^ff — e/cr^ is in- 
dependent of Ac/), while for Hertzian potentials, fee// ~ 

Note that the bulk modulus obeys the same scaling 
0,111, so 



eff ■ 



(28) 



The dimensional length scale in this problem is the 
particle size, a. As a result, the diffusivity is naturally 
measured in units of the product of times the charac- 
teristic frequency scale 



M 



(29) 



Inspection of Fig. [T] reveals that, as oj increases, the dif- 
fusivity decreases rapidly until it saturates at a value 
that we denote by do- N ote f rom Fig. [T] that do ~ 
0.35crV^e///Af = 0.35cr\/e/A/, see Eq. (|29l) . 

Similarly, upon evaluating Eqs. [27] and [21] for a 
Hertzian potential (a = 5/2) in three dimensions, we 
obtain the following prediction for the value of c?o 

d^^c^a^j — Acj?-''^ Hertzian. (30) 



do = c(TyJe/MA4?-^^ where c is of order unity. Fig. 
shows d{uj) for a Hertzian system at four different pack- 
ing fractions. Clearly, do increases with A0. If we divide 
d{u}) and a; by A(j?-/^ as in Fig. [DJb) we get a good col- 
lapse of the data, indicating that do obeys the prediction, 
with c = 0.35. The numerical value of c is independent 
of the potential used and will be derived purely from 
the random geometry of a marginally jammed packing in 
Section [VU 

These results for harmonic and Hertzian potentials 
suggest that the plateau in the diffusivity is consistent 
with what is loosely referred to as the " minimal conduc- 
tivity hypothesis." [13, Hi. The diffusivity do IS mini- 
mal in the sense that the length scale that multiplies the 
characteristic frequency is the smallest length that can be 
chosen in the system, the particle diameter a. Once d{ijj) 
attains this minimal value it cannot decrease by much as 
Lo increases, unless a transition occurs to a new transport 
regime characterized by a vanishingly small diffusivity. 
This is precisely what happens at the end of the plateau 
where Anderson localization sets in. 
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FIG. 2: (Color Online) Plots of diffusivity vs. frequency 
for an unstressed packing of 2000 particles interacting via a 
Hertzian potential. The packing fractions are approximately 
= 0.3 (green), 0.1 (red), 0.05 (blue), 02 (y ellow), (a) 
d{uj) in units of o"y^^) vs. uj (in units of (b) Scaled 

diffusivity, d^uj)/ A<j}^^*, vs. scaled frequency, uij A(^^^^ ^ show- 
ing data collapse in the plateau region. The inset shows the 
packing fraction scaling of the crossover frequency Ud- 



B. Scaling of transport crossover 

The scaling of uj* has been derived for systems near 
the isostatic jamming transition using a variational ar- 
gument that predicts the presence of extended hetero- 
geneous modes with strong spatial de-correlations p^ . 
This structural property suggests that the ability to 
transport heat for vibrational modes above w* should 
be impaired [l6| . In this section, we provide a heuristic 
argument which also suggests that ujd ~ ^* independent 
of the repulsive potential (e.g., for all a). 

Inspection of Fig. [1] and Fig. [2] shows that for the 
range of packing fraction accessible to our simulation, 
the diffusivity appear to be smooth across the frequency 
ujd- At the crossover, where the diffusivity plateau do 
begins and the ballistic regime has just ended, the kinetic 
formula for the diffusivity, Eq. [5] leads 



do 



1 



(31) 



where we have assumed that an effective transverse speed 
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of sound can be assigned to the vibrational modes up to 
the crossover frequency and that it exhibits a smooth 
crossover above it. This working assumption is corrob- 
orated in Section IIVDI where the dispersion relation is 
determined numerically. The phonon density of states at 
low Lo is dominated by shear waves near Point J as appar- 
ent from Equations (jlOjllip and the discussion following 
it. 

We can now solve Eq. PT|) for the mean free path 
id at the crossover and compare it to wavelength 
which we can easily extract from the dispersion relation 
= vtqd- We find that over our range of compression 



(32) 



The crossover frequency, Wrf, corresponds to the "largest" 
wave-vector that can be meaningfully defined in our 
disordered packings and is consistent with the loffe Regel 
criterion [23]. That is to say, the wavelength 1/qd is of 
order of the mean free path id [11] over the range of com- 
pression probed in our study. The main requirement to 
extend the validity of the loffe Regel criterion for smaller 
values of A(/), is that the number on the left-hand side of 
Eq. (j32p remains weakly dependent on packing fraction. 

Upon substituting the mean free path Id from Eq. ([32]) 
into Eq. (|3T|) . we obtain 



1 / ^0 
l/<?d 

Vt 



VKff 



Vt 



(33) 



where the overall scaling of the diffusivity curve indicated 
in Eq. (^5)) was used to write the last step and id ~ l/^d- 
Upon setting the transverse sound speed proportional 
to the square-root of the shear modulus G we find 



Qd 



G 



A(j) 



1/4 



(34) 



where Eqs. Pl9t were substituted for the elastic moduli. 
As ^ — > 0c, the dynamic wavevector qm vanishes with a 
power-law exponent of 1/4, independent of the underly- 
ing inter-molecular potential. This scaling is consistent 
with previous numerical studies of the peak position of 
the transverse structure factor at the frequency onset of 
excess modes, which found a scaling of /\^(j)0-24±0-03 
The non-trivial hypothesis that the diffusivity is smooth 
across tOd implies that id ^ A^^-'^^'*. 

Upon substituting Eq. (j33[) in the dispersion relation 
and using Eqs. [7][S1 we obtain 



G 



VkJT 



Acj)- 



{z - ZcT~ 



(35) 



It is straightforward to conclude from Eq. that LOd 
should scale as A(/)5 and A</)3 for harmonic and hertzian 
interactions respectively. This conclusion is consistent 
with the numerical results plotted in the insets of Figs. [1] 
and[2i;b). 




FIG. 3: Plot of the transport crossover versus the boson peak 
at different Ai^ in jammed packings of 2000 particles with har- 
monic (black open circles) and Hertzian (red open triangles) 
repulsions. The ratio is near unity and is nearly constant with 
A<j!>. 



The limited dynamic range of the simulation stems 
from the fact that the crossover cannot be seen if the 
wavelength l/f/d exceeds the system size. The disappear- 
ance of the low-w diffusivity upturn at low A0 is clearly 
seen in Figs. [TJH In future work, we hope to probe in 
more detail the behavior of the diffusivity in the crossover 
region. 



C. Relation of transport crossover to boson peak 

In this section we show numerically that the transport 
crossover frequency observed in the diffusivity plots 
of Figs. [T][21 corresponds to the same frequency scale, 
w*, at which the onset of the excess vibrational modes is 
observed in the density of states. The signal advantage 
of jammed sphere packings over models studied previ- 
ously is that one can verify this identification at different 
packing fractions and hence test not only that the two 
frequency scales are close in numerical value for a given 
but also that they scale in the same way as a function 
of compression. 

Wc first note that the scaling for Wd in Eq. ([55]) is iden- 
tical with the earlier numerical observation of the boson 
peak frequency, w*, in Eq. (fT2|) as well as with the rela- 
tion derived theoretically in Ref. [IB] for the frequency 
onset of the anomalous modes of compressed jammed 
packings. Fig. [3] shows the ratio ujd/oj* as a function 
of compression. At/). Here, the frequency w* is measured 
numerically from the onset of the plateau in the density 
of states; see Ref. [13] for details. The variation of the 
ratio Ud/oj* is small compared with the variation of w^, 
which changes by an order of magnitude over the same 
range of (see inset to Fig. [Ij. Thus, the transport 
crossover frequency and boson peak frequency track each 
other, implying that the same physics underlies both phe- 
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nomena. In particular, the result shows that the excess 
modes above the boson peak frequency have a small and 
nearly frequency-independent difFusivity. 



D. Change in nature of modes at transport 
crossover 

The Fourier decomposition of the vibrational modes 
evolves dramatically as the frequency is increased 
through the transport crossover at Ud- We concen- 
trate here on fT{q,Lij), the transverse Fourier compo- 
nents [11,11^: 



/T(g,w) = 



^ q A e„(tj) cxp(zq • r„) 



(36) 



where q denotes the wavevector and e„(cj) is the polar- 
ization vector of particle n of the mode at frequency, lo. 
The brackets indicate an average over directions of q. 

At all compressions, fT{q,oj) has a low- wavevector 
peak aX q = qmax that shifts to higher values with in- 
creasing frequency. Two typical examples are shown 
in Figure [4{a) at A0 = 0.1. Below ujd (black curve), 
where the diffusivity decreases rapidly with increasing 
frequency, the peak is sharp and tall. Here the modes 
resemble weakly-scattered transverse plane waves with 
wave- vector qmax- By contrast above uJd (red curve), in 
the region of the diffusivity plateau, the peak is dramat- 
ically less pronounced; the peak height is ~ 50 times 
smaller than for the black curve and is comparable to 
the background signal observed at higher q. This is the 
characteristic signature of the strong-scattering regime 
where vibrational modes can no longer be meaningfully 
characterized by a narrow range of q. Such modes are 
poor conductors of energy, as reflected in the low value 
of d{uj) above cod- This evolution in character is not sharp 
and the peak height decreases continuously as u; increases 
past UJd] a peak - albeit a small one - appears even for 
modes in the diffusivity plateau. 

From data of qmax versus lo one can determine the 
transverse-sound dispersion curve [sl, [5^ . Fig. IDJb-d) 
shows the transverse dispersion relation for three of the 
packing fractions shown in the diffusivity plot of Fig. [TJ 
Each point represents the value of qmax obtained from 
the peak in fT{q,io) for a single vibrational mode of 
frequency to. In other words, each point represents the 
wavevector that makes the largest contribution to a vi- 
brational mode. At each packing fraction the crossover 
frequency ujd is represented by a horizontal dotted line. 
The solid red line, uj = vtq, shows the expected transverse 
dispersion relation with vt decreasing with decreasing 
A0 on the basis of Eq. pH)) . The dashed green line 
in each panel has the same slope independent of Aip. 
From Eq. pO)) . its slope is therefore proportional to (but 
smaller than) the longitudinal sound speed, vi. 

These dispersion curves also show a marked change in 
behavior as the frequency is varied through uJd- Below 



the crossover frequency, the peaks are not only sharp, as 
indicated by the behavior in Figure Ul^a) , but their posi- 
tion corresponds very well with that given by the trans- 
verse speed of sound (red solid line). Above uid, on the 
other hand, the peaks are squat and broad and no longer 
follow the line given by the transverse sound speed. In- 
stead, the lowest values of qmax shift to smaller q and 
begin to track the green line which has a slope indepen- 
dent of compression and proportional to the longitudinal, 
not the transverse, speed of sound. The spread in the po- 
sitions of qmax for frequencies lu > Ud indicate that the 
peaks are very broad in this region. Similar results for 
fT{q,ui) and the transverse dispersion relation above ujd 
were earlier found for a model that included the stress 
terms in the energy [5Q]. In that case, the dispersion re- 
lation above uJd, showed essentially no variation in the 
position of the peaks, qmax, despite an increase in 
by five orders of magnitude. (The data in that study is 
slightly different from the data reported here in that it 
showed averaged /r(9,w) over several nearby frequency 
modes in a bin instead of the peaks in individual modes. 
Both ways of treating the data show the same general 
features.) 

At frequencies below ujd, the transport crossover is her- 
alded by when the weakly scattered plane waves begin to 
scatter strongly. When the modes are only weakly scat- 
tered, the dispersion curves are marked by a single sharp 
peak at q — uj/vt- As the scattering becomes stronger, 
other wave vectors begin to contribute to fxiq, w) so that 
the peaks become very broad and includes a very wide 
range of wave vectors. This is qualitatively what we see 
in Figure Hlja) where the small peak in the red curve is 
hardly greater than the background value. 



V. PLATEAU IN THE DIFFUSIVITY 

In the previous section, we have shown that LUd, where 
the diffusivity flattens out, scales in the same way as does 
the frequency associated with the excess modes of the bo- 
son peak. We can thus consider the flat diffusivity as the 
signature of the excess vibrational modes that appears 
in their transport properties. In the present section, we 
will focus on why the diffusivity is flat over an extended 
frequency range above cud- We will address this (A) by 
showing that from the form of the matrix elements, the 
diffusivity should be simply proportional to the density 
of the modes themselves, which is also flat in this region; 
and (B) by examining what properties of the modes are 
necessary for producing a flat diffusivity. Our study sug- 
gest that ojd ~ UJ* vanishes at Point J. At this point, the 
plateau extends over the entire frequency range up to the 
onset of localization. This simplifies the analysis. Thus, 
Point J is a natural place to gain insight into the origin 
of the constant diffusivity. 
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FIG. 4: (Color online) Transverse mode structure factor 
fT{q,uj) for TV = 2000 and = 0.1 at a; = 0.23 (black) 
below the transport crossover at uim = 0.36 and at = 0.74 
(red) well above the crossover, (b-d) Phono n dispersion re- 
lations for (b) A<;6 = 0.1, (c) 0.05, and (d) 0.01, respectively. 
The data are at discrete values of q because of the finite sys- 
tem size. The horizontal (blue) dotted line in each panel (b- 
d) marks iDd, while the solid line (red) marks the transverse 
sound speed that varies with compression as given by Eq. 1101 
As a comparison, the dashed (green) line is the same in each 
panel (independent of compression) and has a slope propor- 
tional to the speed of longitudinal sound. Note that the solid 
line gives the dispersion for frequencies below tj^. 



A. Energy flux matrix elements 

We start by showing that the energy flux matrix el- 
ements have a particularly simple form at the jamming 
threshold, which enables us to determine the N oo 
behavior of the diffusivity [Tsj. Consider the frequency 
averaged matrix elements defined as 



where i and j are indexes labeling the vibrational modes 
and the matrix elements Sij are given by Eq. (j20p in the 
limit uji LOj. 

Figure [5]^a) shows the heat-flux matrix elements 
|E(tj,a;')P defined in Eq.([37]) for packings at — (/>c = 
10~^ for different values of lo versus lu' . Note that the 
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FIG. 5: (Color Online) Diffusivity just above the jamming 
transition at A(/> = 10^^. (a) Contour plot of the heat- 
flux matrix elements |E(Lj,tj')p plotted versus uj and uj' at 
TV — 2000. (b) Scaling plot showing collapse of |E(tj,aj')p 
at TV = 2000 (black solid), 1000 (red dashed), and 500 (blue 
dotted) with scale factors and w. (c) Scale factors (red 
symbols) and w (black symbols) versus uj. We find oc uP' 
(red dashed line) and w ot lu (black dotted line) except at 
high UJ. id) Scatter Plot of d{u]) with TV = 2000 and delta 
function broadening r] — 0.002 (dashed line). Solid line shows 
predicted d{(jj) for the infinite system. 



matrix elements are symmetric in co and uj' and that they 
increase with increasing uj and uj' . Fig. \^h) shows that 
all the curves for different system sizes, TV, and frequen- 
cies, UJ, can be collapsed onto a simple scaling form given 
by 



\^{UJ,UJ')\' 

|E(c.,c.')P 



1 

TV 
1 

TV 



ui 



if uj' > UJ 
ff u)' <uj (38) 
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except at hi gh f requency where the modes become local- 
ized [HI, [s^nUl and the curves in Fig. E^b) start peel- 
ing off from the green line. Figure EJc) shows that the 
scale factors for the collapse satisfy — uj^ and w — uj, 
respectiyely, consistently with Eq. ((38)) . The scaling col- 
lapse demonstrates that the only noticeable system-size 
dependence is a prefactor of 1/A^. Since for large A^, the 
density of states scales as iV [12], Eq. ^ therefore yields 
a well-defined diffusivity in the N ^ oo limit, shown as 
the solid curve in Fig. [Sfd). 

The scaling collapse in Fig. [IJb) implies that 
oc uj'^/N at low frequencies. This result, com- 
bined with Eq. I19i implies that d(a;) oc D{uj) at low oj. 
This N ^ oo prediction for d{uj) is shown as the solid 
line in Fig. [5^c). Thus, the diffusivity is nearly constant 
down to w = at Point J because the density of states 
is nearly constant there p^ . 

Over most of the frequency range, this N oo pre- 
diction agrees very well with the scattered points in 
Fig.EI^c), which show d{uj) for a system with N=2000. At 
low frequency d{uj) deviates from the solid line and ex- 
hibits an upturn. This upturn is a finite-size artifact that 
arises from replacing the delta- function with a smoothing 
function in Eq. ([2T|) . It scales as uj~^ with a prefactor 
that vanishes as — > oo, 77 — > 0. 



B. Properties of modes in the plateau 

It is important to understand what specific properties 
of the modes give rise to the plateau in the diffusivity. 
To simplify the analysis, we will consider only systems of 
monodisperse particles interacting via harmonic repul- 
sions {a = 2) just above the jamming threshold. The 
argument can be generalized to the bidisperse case stud- 
ied in this paper. 

The starting point for deriving the flat diffusivity and 
for deriving its plateau value, do, is Eq. [^Olfor the matrix 
element of the heat flux operator evaluated using periodic 
boundary conditions. Recall that for two modes of fre- 
quencies within a bin centered at lu the matrix element 
Sij reads 

Sij = ■^—^e{n,i)H{x,y)e{m,j){Rjn-Rn) (39) 

where n and m label the particles, i?„ and R„i their po- 
sitions, H(n,m) is the dynamical matrix element (itself 
a, d * d matrix) between these two particles, e(n, i) is the 
displacement of particle n in mode i and M is the particle 
mass. 

Next set (Rm — Rn) — o'Rnm, where a is the particle 
diameter and Rnm the unit vector from n to m. Set the 
non-diagonal terms H{n, m) — kRnm Rnm, where k is 
the contact stiffness. With this substitution, Eq. (IHS)) 



can be rewritten as a sum on all contacts (n, m) : 
Sij = X! Rn7n[{e{n,i) ■ Rnm){e{m,j) ■ Rn,n) 

{n,rn) 

- {e{n,j) ■ Rnm){e{m,i) ■ Rnm)] (40) 
We can then rewrite the sum in Eq. (PD|) as: 

T,{uj) = ^ Rnrn[iein,i) ■ R„m - e{m,i) ■ Rnrn)ieirnJ) ■ R„rn) 

{n,m) 

- (e(n, j) ■ Rnrn " e(r7i, j) ■ Rnm)ie{m, i) ■ Rnm)] (41) 

and 

{n,m) 

- {5ri^){e{m,i) ■ Rnm)] (42) 

where <5r^„ = (e(n, i) ■ Rnm — e(m, i) ■ Rnm) is the stretch- 
ing of the contact nm corresponding to mode i. Taking 
the amplitude squared of E leads to diagonal and non- 
diagonal terms. The latter are of two forms: 

{SrnmRnm)(Sr'pgRpq){e{m,j) ■ Rnm){e{q,j) ■ Rpq) (43) 

{^rnmRnm){SrlqRpq){e{m,i) ■ ^„„)(e(p,z) • Rpq) (44) 

where nm and pq correspond to two distinct contacts. 

We now make some assumptions about the nature of 
the modes whose validity we test numerically, (i) if the 
displacements are uncorrelated between different modes 
((e(m,j) • e{q,i)) — 0), the two terms above become: 

{iSrnmRnm){Srlqnpg)){{e{m,j) ■ Rnm){eiq,j) ■ i?pg)l45) 
{{Srl^mRnm){e{q,i) ■ i?p,)) ((e(m, j) • i?„„,)(fr^,i?p,)l46) 

This assumption was tested numerically in a packing 
comprised of iV = 2000 particles at A(f> = 10~^ for the 
low frequency modes. Each of the two terms was found 
to be of the order of 10~^^ which is vanishingly small 
within numerical precision. 

We next assume that modes have weak spatial cor- 
relations. Numerical simulations actually support that 
such correlations exist [s^ and grow as the frequency de- 
creases. Neglecting them nevertheless appears to yield 
good quantitative results, as we shall see below. In par- 
ticular, each of these two terms in Eq. lj45p vanish un- 
der the specific assumptions that (ii) the directions of 
the stretching of different contacts within a mode are 
not correlated in space {{{5rnniRnm){5rpqRpq)) = 0) and 
that (iii) the directions of stretching and of displacement 
are locally uncorrelated (((5r^„^„m)(e(g, i) ■ Rpq)) — 0). 
For assumption (ii) we find that the corresponding terms 
are of the order of 10^^^ or lower. For assumption (iii) 
we find a term of the order of 10~^^ or less. If the quan- 
tities were correlated, we should find that for extended 
modes they are of the order of 1/N k, 10~*, which is 
much larger than the values shown above. Therefore, the 
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assumptions listed above appear to be reasonable and we 
are left with: 

« 2 SrU{e{m,j) ■ R^^f (47) 

(n,m) 

Assuming now that (iv) the amplitude of the displace- 
ments between two modes are not correlated (which is 
wrong for localized modes, where the displacements are 
anti-correlated since localized modes do not live on the 
same regions of space) then we have: 

|Sp = 2 ^ (K™')((e(m, j) . i?„„0^) (48) 

(n,m) 

The modes are normalized so {{e{m, j) ■ Rnm)'^) — l/NV, 
where V is the dimensionality of space. Also, the mode 
energy is SE 
tain 



M^V2 ^ fc/2E(„,™)K,n', so we ob- 



2Muj^ 1 



NV 



Then we have for Eq. (|iD|) : 



This leads to: 



d{uj) 



2V^VN 



TrD{Lu)a^k 
3mV 



(49) 



(50) 



(51) 



where D{ll!) is the density of states per particle. In units 
where cr = fc = m = l, we obtain 



do = 



(52) 



in three dimensions, where dp and Dq denote the plateau 
values of the diffusivity and the density of states, respec- 
tively. This is consistent with the numerical data in Fig. 
[5jd) , which shows that the diffusivity roughly follows the 
nearly flat density of states at Point J. 

In summary, we obtain a frequency-independent diffu- 
sivity when the density of states is frequency-independent 
and the following conditions are satisfied: (A) displace- 
ments of particles in different modes of similar frequency 
are uncorrelated; (B) the directions of changes in the rel- 
ative displacements of pairs of interacting particles are 
spatially uncorrelated within a given mode; and (C) the 
direction of change in the relative displacement of a pair 
of interacting particles is uncorrelated from the direction 
of the displacement. 



C. Stressed packings 

Until now, we have neglected the forces that particles 
within a jammed packing exert on each other by replac- 
ing compressed springs between particles with springs at 




FIG. 6: (Color Online) Diffusivity for A(fi = 0.1 in stressed 
(black) and unstressed (red) packings composed of 1000 par- 
ticles with harmonic repulsions. 



their equilibrium length (see Sec. IIIIBp . Here we restore 
the stress terms into the dynamical matrix and examine 
how they affect the behavior. 

The repulsive forces between particles tend to desta- 
bilize the system towards buckling, or rearrangements. 
As a result, they tend to push modes to lower oj. As a 
result, finite-size effects, which cut off plane waves at low 
frequency, are more obstructive. This prevents us from 
studying the transport crossover directly in stressed sys- 
tems. 

Nevertheless, we suggest that the scaling of ujd with 
compression should be the same in the stressed case as 
in the unstressed one. This follows from the fact the 
boson peak frequency, to*, follows the same power law 
(Eq. [121) as in the unstressed case. This is illustrated in 
Fig. [3 which shows that the low- frequency portion of the 
vibrational spectrum collapses onto a single curve when 
u) is scaled by uj* ^ A0^/^ for harmonic repulsions. 



VI. AC THERMAL CONDUCTIVITY AT 
POINT J. 



The calculation of the energy flux matrix elements in 
Eq. ([55)1 opens up the possibihty of estimating the ther- 
mal conductivity k(T, fl) of the marginally jammed solid 
in the presence of an AC thermal gradient driven with 
angular frequency fi. As a starting point we rewrite Eq. 
(|16p as an integral over uj rather than a double sum over 
eigenmodes 



K{T,n) = - 



nV 



Tfl 

X D{uj) D{uj + n) duj 



n{uj + VL)- n{uj)] \ {lo + n\S\uj)\^ x 

(53) 
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FIG. 7: (Color Online) Density of states D(iJ) vs. oj/oj* for 
Acj!> = 0.0001 (black dotted), 0.01 (yellow dashed-dot-dot), 
0.05 (blue dot-dashed), 0.1 (red dashed) and 0.3 (green solid), 
for a stressed system with = 2000 with harmonic repul- 
sions. Note that tj* scales according to Eq. 1121 This scaling 
produces a good collapse of the vibrational spectrum at low 
frequency. 



where the frequency averaged matrix elements read 

(2w + 0)2 



4tj(cj + ?L) 

uP- ( 2Vl\ 
AT 1 + - 



|E(17 + cj,w)|^ (54) 



(55) 



and higher order corrections in terms of fi^ have been 
dropped. The prefactor in Eq. ([20|) . which was set to 
unity in evaluating S(Li;,a;'), has been explicitly included 
since the mode coupling between w and w' is no longer 
restricted to vibrational states at the same frequency. 
Despite the concise notation adopted, Eq. ([55)) accounts 
for both upwards and downwards jumps, w ^ w ± SI , 
corresponding to energy being absorbed or injected into 
the reservoirs [60|. In order to obtain Eq. (|55p. Eq. ([55]) 
was substituted into Eq. (f54|) . 

With the aid of Eq. (fT5|) , we can expand the difference 
in occupation numbers to first order in O 



T 

hijj' 
d_ 



T 



(56) 



n{u! + fl) — n{u!) 



Upon substituting Equations ([55)1 and ()55p into Eq. ((55 
four terms are obtained of which one is 0{il^) and it will 
be ignored. After performing an integration by parts on 
the 0(0) term involving the lo derivative and canceling 
out two terms which are equal and opposite, we obtain 



h k{T, n) 
ttcqN^ 



C(W) dLU + n [C{lu)]q' 



(57) 



where Nq is the value of the plateau in the density of 
states. 



The desired result follows, according to Eq. 
setting C{uJmax = oo) « 



([T8| . upon 



k{t, n) 



a T — - — 

Kb 



+ 0(^2) , (58) 



where the numerical constant a is given by 



dx 



(e^ - 1)' 



(59) 



fre- 



As expected intuitively, a non-vanishing driving 
quency O results in a lower thermal conductivity. 

Eq. [55] allows the calculation of the T-dependent DC 
thermal conductivity (O 0) at point J . We note that 
this is particularly simple and can be calculated at the 
harmonic level without facing any divergences because 
no acoustic phonons are present. From Eq. 1581 we find 
that k[T) grows linearly in temperature for small T and 



saturates above a temperature kBT„. 



fku„ 



where 



^max is the maximum frequency above which there are 
no more vibrational states. This follows from assuming 
that both the density of states and the diffusivity are 
approximately uj independent as expected from the scal- 
ing analysis and numerical extrapolations presented in 
Section flV Al Thus, a plateau in the diffusivity leads to 
an approximately linear increase of k{T) followed by a 
saturation. 



VII. 



CONCLUSION 



We have studied a class of model amorphous solids 
whose elastic properties can be tuned by varying the den- 
sity near the jamming/ unjamming transition. Point J. 
The proximity to Point J allows variation of the crossover 
frequency that marks the onset of the plateau in the dif- 
fusivity, ujd, with A(j). As A(j) — )■ 0, our scaling arguments 
show that ujd — *■ 0, so that the plateau in the diffusivity 
extends all the way down to zero frequency. Moreover, 
the value of ujd agrees with the boson peak frequency w* 
within a factor close to unity; they both scale the same 
way with A0. 

The findings presented in this study enable us to es- 
tablish that (1) there is a frequency regime in which the 
diffusivity is small and nearly constant, and (2) the bo- 
son peak frequency coincides with the energy transport 
crossover frequency for all pressures applied to our un- 
stressed amorphous packings of repulsive spheres. More 
work is needed to assess the relationship between the bo- 
son peak and the transport cross-over when pre-stress is 
important. Nonetheless, our results suggest solutions to 
two longstanding conundrums posed by these amorphous 
solids. First, in such systems, the plateau in the thermal 
conductivity at intermediate temperatures is followed by 
a rise and then a saturation at high temperatures. Crys- 
tals show the opposite behavior, with a thermal conduc- 
tivity that decreases with T at high T, see Fig. 5.1 in 
Ref. What is the origin of the generic monotonicity 
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of K with T in amorphous solids? Second, the tempera- 
ture range of the plateau in the thermal conductivity lies 
near the temperature at which the heat capacity exhibits 
a boson peak. Is this a coincidence? 

The answer to the first conundrum follows naturally 
from our result that there is a frequency regime of small 
and constant diffusivity. As shown in Sec. lVIi the plateau 
in the diffusivity above LOd gives rise to a linear increase 
in thermal conductivity above ksTd ~ hjJd- Similarly, 
the vanishing of the density of states at oJmax leads to 
saturation of k{T) at high fc^T > hiOmax- Thus, the 
transport crossover frequency sets the high-temperature 
limit of the plateau in k(T), while the maximum allowed 
frequency sets the temperature at which k(T) saturates 
to its final high-temperature value. 

The answer to the second conundrum follows directly 
from the observation that ujd ~ 'j-'*, since the high- 
temperature limit of the plateau in k(T) is ksTd ^ tuod 
and the boson peak temperature is ksTsp ^ fko*. 
Therefore, Td~TBp. 

One might object that our results were obtained for a 
special system in which spheres interact via finite-ranged 
repulsions. It has been argued theoretically [3] and 
shown numerically js^ that packings of spheres interact- 
ing via Lennard- Jones interactions, which are attractive 
at long distances, behave much as compressed packings 
of repulsive spheres and that the boson peak frequency 
shifts upwards with increasing density in such systems. 
Indeed, it was found that the lo* is determined primarily 
by the repulsive interactions that come into play because 
the systems are held at high densities by the attractions. 
Thus, even systems with attractions display a boson peak 
corresponding to the onset of anomalous modes, which 
should lead to a transport crossover as well. 

It has also been shown that packings of ellipsoids [6lj 
display boson peaks corresponding to the onset of anoma- 



lous modes that are similar in character to those for 
spheres. For ellipsoids, the boson peak frequency is con- 
trolled by much the same physics as for spheres; it de- 
pends on the coordination number as in Eq. psp . 

Another class of systems, network glasses, appears at 
first glance to be very different from our model repul- 
sive sphere packings. However, the unstressed models 
that were the main focus of this paper can be described 
as points (corresponding to the centers of spheres) con- 
nected by unstretched springs. Such networks bear 
some resemblance to network glasses, which are held to- 
gether by covalent attractions. We include only central 
forces, but the counting of constraints has been shown 
to be key to network glasses with bond-bending forces as 
well [6^ [63| . The connection between harmonic spring 
networks and covalent network glasses has been discussed 
in more detail in Refs. [l^ and [63 |. 

In summary, our results suggest that two generic fea- 
tures of amorphous solids, the rise of the thermal con- 
ductivity with temperature above a plateau and the co- 
incidence of the plateau temperature with the boson 
peak temperature, can be viewed as echoes of Point J. 
This transition controls the onset frequency of anomalous 
modes with constant, minimal diffusivity in the manner 
of a critical point. 
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